Objetivo: trabajar con algoritmos de machine learning en R para la predicción espacio-temporal de variables.
Machine learning es una rama de la Inteligencia Artificial que tiene como objetivo utilizar datos y algoritmos para imitar el aprendizaje de los seres humanos (IBM), y se enfoca en realizartareas conocidas.
Los algoritmos de machine learning se han utilizado ampliamente en:
El machine learning ha demostrado ser útil en muchos campos. ¿Cuáles son los pasos para aplicar métodos de machine learning?
No existe un algoritmo de machine learning que sea el mejor en todos los casos. Su rendimiento depende del tipo de datos y los objetivos finales. Básicamente, los dos principales grupos de algoritmos de machine learning son:
Los siguientes algoritmos se utilizan a menudo para el aprendizaje supervisado:
Los siguientes algoritmos se utilizan a menudo para el aprendizaje no supervisado:
En este curso, nos centraremos principalmente en algoritmos de aprendizaje supervisado, ya que queremos trabajar con regresión y clasificación de variables.
Random Forest (RF; Biau and Scornet, 2016, Breiman, 2001, Prasad et al., 2006) es un método que se puede utilizar para tareas de clasificación y regresión supervisadas construyendo árboles de decisión utilizando la relación entre variables independientes y dependientes.
\[ \hat{\theta}^B(x) = \frac{1}{B} \sum_{b=1}^B{t^*_b(x)} \] \(\hat{\theta}^B(x)\) es la predicción final
\(b\) es la muestra individual de arranque
\(B\) es el número total de árboles
\(t^*_b\) es el árbol de decisión individual
Support Vector Machines (SVM; Hearst et al., 1998, Steinwart et al., 2008) tiene como objetivo encontrar un hiperplano en un espacio de N dimensiones que clasifique distintamente los puntos de datos dados. Como hay muchos hiperplanos que se pueden seleccionar, SVM maximiza la distancia de margen entre puntos de datos de diferentes clases.
Las Redes Neuronales (NN), también conocidas como Redes Neuronales Artificiales, están compuestas por capas de nodos. Tienen una capa de entrada, una o más capas ocultas y una capa de salida. Cada nodo se conecta a otro y tiene un peso y un umbral asociados (IBM). Si la salida de un nodo diseñado está por encima del umbral especificado, el nodo envía información al siguiente.
Ahora que hemos obtenido una visión general de diferentes algoritmos de machine learning, vamos a usarlos en un ejercicio de regresión. El conjunto de datos mtcars (incluido en R) comprende el consumo de combustible y 10 aspectos del diseño y rendimiento de automóviles para 32 automóviles.
print(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2Queremos predecir la potencia (hp) de algunos autos
utilizando los otros parámetros de los autos:
En este ejercicio utilizaremos Random Forest, Support Vector Machines y Neural Networks. Por lo tanto, carguemos los paquetes correspondientes:
Para predecir la potencia (hp) de los autos, vamos a separar aleatoriamente los datos en una muestra de entrenamiento (para entrenar los algoritmos de machine learning) y una muestra de validación (para evaluar su rendimiento). Para garantizar la reproducibilidad, estableceremos una semilla en el código. Establecer una semilla en R significa inicializar un generador de números pseudoaleatorios.
Para separar los datos, vamos a almacenar en un objeto el número total de autos.
Después de esto, podemos generar una muestra aleatoria de posiciones. Vamos a utilizar el 80% del número total de autos con fines de entrenamiento.
Ahora podemos separar el conjunto de datos mtcars en conjuntos de entrenamiento y validación.
training <- mtcars[pos_training,]
validation <- mtcars[-pos_training,]
print(training)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
print(validation)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1El paquete randomForest implementa el algoritmo de random forest de Breiman para clasificación y regresión.
x: data frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar
data: data frame que contiene las variables del modelo
ntree: número de árboles a construir
mtry: número de variables muestreadas aleatoriamente como candidatos en cada división
nodesize: Tamaño mínimo de los nodos terminales. Si se establece un número mayor, se generarán árboles más pequeños (y, por lo tanto, se tardará menos tiempo)
maxnodes: Número máximo de nodos terminales que pueden tener los árboles del bosque. Si no se especifica, los árboles se crecerán hasta el máximo posible (sujeto a los límites de nodesize)
na.action: Una función para especificar la acción a tomar si se encuentran valores faltantes (NA)
Existen más parámetros relacionados con esta función. Por favor, escribe ?randomForest para verlos.
En este ejemplo, utilizaremos los valores predeterminados de la función. La fórmula que se aplicará es \(hp\sim.\). Esto significa que \(hp\) es la variable dependiente y el \(.\) significa que se utilizarán todas las variables del conjunto de datos de entrenamiento.
En caso de que queramos utilizar solo parámetros específicos, podemos configurarlos de la siguiente manera:
Podemos imprimir el objeto rf_model.
print(rf_model)
##
## Call:
## randomForest(formula = hp ~ ., data = training)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 1322.719
## % Var explained: 73.54Para predecir los valores de hp, utilizaremos la función predict, que requiere el modelo y las covariables para predecir la potencia de los automóviles. Para dejarlo explícito, excluiremos la columna de hp del data.frame.
Podemos crear un data.frame con los valores reales (Obs) y las predicciones del Random Forest (RF).
La función varImpPlot generará un gráfico de puntos de la importancia de las variables según lo medido por un Random Forest.
El paquete e1071 se utiliza para entrenar una máquina de vectores de soporte para clasificación y regresión.
formula: un data.frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar
data: un data.frame que contiene las variables del modelo
kernel: el núcleo utilizado en el entrenamiento y la predicción: lineal, polinómico, radial y sigmoidal
degree: parámetro relacionado con el núcleo polynomial
gamma: parámetro necesario para todos los núcleos excepto el lineal
coef0: parámetro necesario para los núcleos de tipo polinómico y sigmoidal
na.action: una función para especificar la acción a tomar si se encuentran valores faltantes (NA)
Existen más parámetros relacionados con esta función. Por favor, escribe ?svm para verlos.
Al igual que RF, la fórmula que se aplicará es \(hp\sim.\). Utilizaremos diferentes núcleos.
svm_model_lin <- svm(hp ~ ., data = training, kernel = "linear")
svm_model_poly <- svm(hp ~ ., data = training, kernel = "polynomial")
svm_model_rad <- svm(hp ~ ., data = training, kernel = "radial")
svm_model_sig <- svm(hp ~ ., data = training, kernel = "sigmoid")Ahora podemos predecir los valores de hp:
Finalmente, podemos resumir los resultados.
data.frame(Obs = validation[,4], SVM_lin = svm_prediction_lin,
SVM_poly = svm_prediction_poly, SVM_rad = svm_prediction_rad,
SVM_sig = svm_prediction_sig)
## Obs SVM_lin SVM_poly SVM_rad SVM_sig
## Hornet 4 Drive 110 121.51574 123.47778 106.47474 120.99467
## Merc 280 123 156.98698 146.04602 125.07476 142.84838
## Merc 450SL 180 156.65475 168.71280 171.96696 160.63893
## Cadillac Fleetwood 205 236.92939 223.57702 211.27667 226.52371
## Toyota Corona 97 45.72309 80.75757 91.49229 82.29314
## Fiat X1-9 66 69.46686 77.96054 77.80608 68.92851El paquete neuralnet permite entrenar redes neuronales utilizando retropropagación, retropropagación resiliente (RPROP) con o sin retroceso de peso.
formula: un data.frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar
data: un data.frame que contiene las variables del modelo
hidden: un vector de enteros que especifica el número de neuronas ocultas en cada capa
threshold: un valor numérico que especifica el umbral para las derivadas parciales de la función de error como criterio de detención
stepmax: el número máximo de pasos para el entrenamiento de la red neuronal
algorithm: una cadena de texto que contiene el tipo de algoritmo para calcular la red neuronal. Los tipos posibles son: ‘backprop’, ‘rprop+’, ‘rprop-’, ‘sag’ o ‘slr’
Existen más parámetros relacionados con esta función. Por favor, escribe ?neuralnet para verlos.
La fórmula que se aplicará es \(hp\sim.\).
Ahora podemos predecir los valores de hp:
Finalmente, podemos resumir los resultados.
data.frame(Obs=validation[,4], NN=nn_prediction)
## Obs NN
## Hornet 4 Drive 110 150.4998
## Merc 280 123 150.4998
## Merc 450SL 180 150.4998
## Cadillac Fleetwood 205 150.4998
## Toyota Corona 97 150.4998
## Fiat X1-9 66 150.4998La primera razón a considerar cuando se obtienen resultados extraños con redes neuronales es la normalización. Los datos deben estar normalizados, de lo contrario, sí, el entrenamiento dará como resultado una NN sesgada que producirá el mismo resultado todo el tiempo.
Normalicemos los datos.
norm_df <-rbind(training, validation)
#define Min-Max normalization function
minmax_norm <- function(x) {
(x - min(x)) / (max(x) - min(x))
}
norm <- apply(norm_df, 2, minmax_norm)
head(norm)
## mpg cyl disp hp drat wt
## Chrysler Imperial 0.1829787 1.0 0.9201796 0.6289753 0.2165899 0.9798006
## Volvo 142E 0.4680851 0.0 0.1244699 0.2014134 0.6221198 0.3239581
## Merc 230 0.5276596 0.0 0.1738588 0.1519435 0.5345622 0.4185630
## Maserati Bora 0.1957447 1.0 0.5734597 1.0000000 0.3594470 0.5259524
## Mazda RX4 Wag 0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.3482485
## Merc 450SLC 0.2042553 1.0 0.5106011 0.4522968 0.1428571 0.5796471
## qsec vs am gear carb
## Chrysler Imperial 0.34761905 0 0 0.0 0.4285714
## Volvo 142E 0.48809524 1 1 0.5 0.1428571
## Merc 230 1.00000000 1 0 0.5 0.1428571
## Maserati Bora 0.01190476 0 1 1.0 1.0000000
## Mazda RX4 Wag 0.30000000 0 1 0.5 0.4285714
## Merc 450SLC 0.41666667 0 0 0.0 0.2857143Ahora, podemos seleccionar los datos correspondientes.
Finalmente, podemos aplicar el modelo de NN.
Podemos graficar la NN.
Podemos graficar la NN de una manera más agradable.
Podemos usar los modelos para predecir la variable.
nn_prediction_log <- predict(nn_model_log, norm_validation[,-4])
nn_prediction_tan <- predict(nn_model_tan, norm_validation[,-4])Para visualizar los resultados:
data.frame(Obs = norm_validation[,4], NN_log = nn_prediction_log, NN_tan = nn_prediction_tan)
## Obs NN_log NN_tan
## Hornet 4 Drive 0.20494700 0.15848709 0.280470546
## Merc 280 0.25088339 0.30682250 0.338217726
## Merc 450SL 0.45229682 0.35768535 0.359986470
## Cadillac Fleetwood 0.54063604 0.61909455 0.705239030
## Toyota Corona 0.15901060 0.10892604 -0.003904399
## Fiat X1-9 0.04946996 0.08748578 0.041047591Finalmente, podemos volver a transformar los datos de hp:
hp_val_log <- nn_prediction_log * (max(norm_df$hp) - min(norm_df$hp)) + min(norm_df$hp)
hp_val_tan <- nn_prediction_tan * (max(norm_df$hp) - min(norm_df$hp)) + min(norm_df$hp)
data.frame(Obs = validation[,4], NN_log = hp_val_log, NN_tan = hp_val_tan)
## Obs NN_log NN_tan
## Hornet 4 Drive 110 96.85185 131.37316
## Merc 280 123 138.83077 147.71562
## Merc 450SL 180 153.22495 153.87617
## Cadillac Fleetwood 205 227.20376 251.58265
## Toyota Corona 97 82.82607 50.89506
## Fiat X1-9 66 76.75848 63.61647Finalmente, creemos un resumen de los resultados.
result <- data.frame(Obs = validation[,4], RF = rf_prediction, SVM_lin = svm_prediction_lin,
SVM_poly = svm_prediction_poly, SVM_sig = svm_prediction_sig,
NN_log = hp_val_log, NN_tan = hp_val_tan)
print(result)
## Obs RF SVM_lin SVM_poly SVM_sig NN_log
## Hornet 4 Drive 110 123.43118 121.51574 123.47778 120.99467 96.85185
## Merc 280 123 122.20713 156.98698 146.04602 142.84838 138.83077
## Merc 450SL 180 173.63778 156.65475 168.71280 160.63893 153.22495
## Cadillac Fleetwood 205 219.34986 236.92939 223.57702 226.52371 227.20376
## Toyota Corona 97 96.53314 45.72309 80.75757 82.29314 82.82607
## Fiat X1-9 66 74.71494 69.46686 77.96054 68.92851 76.75848
## NN_tan
## Hornet 4 Drive 131.37316
## Merc 280 147.71562
## Merc 450SL 153.87617
## Cadillac Fleetwood 251.58265
## Toyota Corona 50.89506
## Fiat X1-9 63.61647Ahora hagamos un ejemplo más práctico. Imagina que tenemos dos cuencas anidadas y tenemos algunos datos faltantes en la cuenca inferior que necesitamos. Usemos el aprendizaje automático para predecir esos valores faltantes.
Para trabajar con la serie temporal, utilizaremos el paquete hydroTSM. Primero, vamos a leer los datos, que se almacenan como un objeto zoo.
Separemos los datos de la cuenca y la subcuenca en consecuencia.
Dado que en este caso tenemos los datos completos, crearemos un objeto adicional con el caudal de la cuenca y eliminaremos algunos años de datos. De esta manera, podremos comparar los valores medidos con las predicciones.
Podemos almacenar las fechas del objeto zoo en un objeto.
Ahora, estableceremos los primeros seis años de catch_na como NAs. Para este propósito, crearemos un objeto pos que almacenará la posición de los datos relacionados con ‘1995-12-31’.
En este sentido, el objeto catch tendrá los datos completos.
Mientras que el objeto catch_na tendrá seis años de valores faltantes.
Particionaremos los datos que se utilizarán para fines de entrenamiento y verificación utilizando el objeto pos que creamos anteriormente.
# El conjunto de entrenamiento debe comenzar un día después de "1995-12-31"
inicio <- pos + 1
subcatch_training <- subcatch[, inicio:length(subcatch)]
catch_training <- catch[, inicio:length(catch)]
head(subcatch_training)
## 1990-01-01 1990-01-02 1990-01-03 1990-01-04 1990-01-05 1990-01-06
## 4.540620 4.428097 4.269241 4.203052 4.090529 4.004482
head(catch_training)
## 1990-01-01 1990-01-02 1990-01-03 1990-01-04 1990-01-05 1990-01-06
## 2.521102 2.452117 2.392539 2.389403 2.285925 2.220075training <- data.frame(subcatch = subcatch_training, catch = catch_training)
head(training)
## subcatch catch
## 1990-01-01 4.540620 2.521102
## 1990-01-02 4.428097 2.452117
## 1990-01-03 4.269241 2.392539
## 1990-01-04 4.203052 2.389403
## 1990-01-05 4.090529 2.285925
## 1990-01-06 4.004482 2.220075
tail(training)
## subcatch catch
## 2018-12-26 3.263157 1.859470
## 2018-12-27 3.203586 1.803027
## 2018-12-28 3.130777 1.774806
## 2018-12-29 3.177110 1.771670
## 2018-12-30 3.276395 1.909641
## 2018-12-31 3.005016 1.762263
summary(training)
## subcatch catch
## Min. : 0.9664 Min. : 0.4829
## 1st Qu.: 2.7403 1st Qu.: 1.3672
## Median : 4.5208 Median : 2.8190
## Mean : 5.7265 Mean : 4.1321
## 3rd Qu.: 7.2147 3rd Qu.: 5.7070
## Max. :45.2076 Max. :47.0982
## NA's :123 NA's :7Durante el entrenamiento, podemos enfrentar algunos problemas si tenemos datos faltantes. Excluyamos los datos faltantes del entrenamiento.
pos_nas <- which(is.na(training$subcatch))
pos_nas <- c(pos_nas, which(is.na(training$catch)))
print(pos_nas)
## [1] 7396 7397 7398 7399 7400 7401 7402 7403 7404 7405 7406 7407 7408 7409 7410
## [16] 7411 7412 7413 7414 7415 7416 7417 7418 7419 7420 7421 7422 7423 7424 7425
## [31] 7426 7427 7428 7429 7430 7431 7432 7433 7434 7435 7436 7437 7438 7439 7440
## [46] 7441 7442 7443 7444 7445 7446 7447 7448 7449 7450 7451 7452 7453 7454 7455
## [61] 7456 8188 8189 9442 9443 9444 9445 9446 9447 9448 9449 9450 9451 9453 9454
## [76] 9455 9456 9457 9458 9460 9461 9462 9463 9464 9465 9466 9467 9468 9469 9470
## [91] 9471 9472 9473 9474 9475 9476 9477 9478 9479 9480 9481 9482 9483 9484 9485
## [106] 9486 9487 9488 9489 9490 9491 9492 9493 9494 9495 9496 9497 9498 9499 9798
## [121] 9799 9800 9801 2237 2238 2239 2240 5537 6686 9831
training <- training[-pos_nas,]
summary(training)
## subcatch catch
## Min. : 0.9664 Min. : 0.4829
## 1st Qu.: 2.7403 1st Qu.: 1.3609
## Median : 4.5274 Median : 2.8692
## Mean : 5.7287 Mean : 4.1585
## 3rd Qu.: 7.2147 3rd Qu.: 5.7383
## Max. :45.2076 Max. :47.0982El siguiente paso es entrenar los modelos. Para este propósito, utilizaremos Random Forest y Support Vector Machines con los diferentes kernels.
model_rf <- randomForest(catch ~ ., data = training)
model_svm_lin <- svm(catch ~ ., data = training, kernel = "linear")
model_svm_poly <- svm(catch ~ ., data = training, kernel = "polynomial")
model_svm_rad <- svm(catch ~ ., data = training, kernel = "radial")
model_svm_sig <- svm(catch ~ ., data = training, kernel = "sigmoid")Ahora, preparemos los datos que se utilizarán para predecir los valores para la subcuenca (es decir, el caudal diario de los primeros seis años de la cuenca).
covariate <- data.frame(subcatch = subcatch[1:pos])
head(covariate)
## subcatch
## 1990-01-01 4.540620
## 1990-01-02 4.428097
## 1990-01-03 4.269241
## 1990-01-04 4.203052
## 1990-01-05 4.090529
## 1990-01-06 4.004482
tail(covariate)
## subcatch
## 1995-12-26 3.918436
## 1995-12-27 3.819151
## 1995-12-28 3.653676
## 1995-12-29 3.541154
## 1995-12-30 3.441869
## 1995-12-31 3.375679Finalmente, podemos predecir los valores de caudal utilizando los diferentes modelos y el objeto covariable.
Ahora que tenemos los valores predichos, podemos convertirlos en objetos zoo.
Podemos graficar los resultados de las predicciones de la siguiente manera:
plot(catch_na)
lines(prediction_rf, col="green")
lines(prediction_svm_lin, col="blue")
lines(prediction_svm_poly, col="red")
lines(prediction_svm_rad, col="cyan")
lines(prediction_svm_sig, col="purple")
grid()
legend("topright", lwd = 1, c("RF", "SVM-Lineal", "SVM-Polinomial", "SVM-Radial", "SVM-Sigmoid"),
col = c("black", "green", "blue", "red", "cyan", "purple"), bty = "n")Podemos graficar los resultados de las predicciones de la siguiente
manera:
Eliminando los modelos de bajo rendimiento:
Podemos graficar los resultados de las predicciones de la siguiente
manera:
Podemos utilizar un indicador de rendimiento para evaluar qué modelo tuvo el mejor desempeño. La eficiencia de Kling-Gupta (KGE; Gupta et al., 2009, Kling et al., 2012) tiene tres componentes: el coeficiente de correlación de Pearson (\(r\)); la relación de sesgo (\(\beta\)); y la relación de variabilidad (\(\gamma\)).
\[\scriptsize \textrm{KGE'} = 1 - \sqrt{ (r-1)^2 + (\beta-1)^2 + (\gamma-1)^2 } \]
\[\scriptsize r = \frac{ \sum_{i=1}^n{(O_i - \bar{O}) (S_i - \bar{S}) } }{\sqrt{\sum_{i=1}^n{(O_i - \bar{O})^2}}\sqrt{\sum_{i=1}^n{(S_i - \bar{S})^2}}}; \beta = \frac{\mu_{s}}{\mu_{o}}; \gamma = \frac{CV_{s}}{CV_{o}} = \frac{\sigma_{s}/\mu_{s}}{\sigma_{o}/\mu_{o}} \]
Donde \(\mu\) es el promedio del caudal (\(Q\)), \(CV\) es el coeficiente de variación, \(\sigma\) representa la desviación estándar de \(Q\), y los subíndices \(s\) y \(o\) representan el caudal simulado y observado, respectivamente. El KGE’ y sus componentes tienen su valor óptimo en uno. \(r\) se relaciona a la dinámica temporal, mientras \(\beta\) y \(\gamma\) se relacionan a el volumen y la variabilidad, respectivamente.
Podemos crear un objeto que almacene las observaciones del período de verificación.
Finalmente, podemos aplicar el KGE’ utilizando la función KGE del paquete hydroGOF.
KGE(obs = obs, sim = prediction_rf, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.9234464
##
## $KGE.elements
## r Beta Gamma
## 0.9682697 1.0035594 0.9304229
KGE(obs = obs, sim = prediction_svm_lin, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.9080331
##
## $KGE.elements
## r Beta Gamma
## 0.9559886 0.9970294 0.9193025
KGE(obs = obs, sim = prediction_svm_poly, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.08816593
##
## $KGE.elements
## r Beta Gamma
## 0.6845051 0.7331579 1.8128343
KGE(obs = obs, sim = prediction_svm_rad, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.8847093
##
## $KGE.elements
## r Beta Gamma
## 0.9517110 0.9925498 0.8955749
KGE(obs = obs, sim = prediction_svm_sig, method = "2012", out.type = "full")
## $KGE.value
## [1] -90.43503
##
## $KGE.elements
## r Beta Gamma
## -0.7320031 -90.2282000 -4.8975169